\documentclass{article}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{geometry}
\usepackage{hyperref}
\usepackage{graphicx}
\usepackage{float}
\usepackage{subcaption}

\graphicspath{{figures/}}

\geometry{a4paper, margin=1in}

\title{A Numerical Framework for Simulating the Aerodynamics of Dragonfly Flight}
\author{Your Name \\ Your Institution}
\date{\today}

\begin{document}

\maketitle 

\begin{abstract}
The remarkable flight capabilities of dragonflies, characterized by high lift, exceptional maneuverability, and hovering, have long inspired the design of Micro Air Vehicles (MAVs). Understanding the complex, unsteady aerodynamics underlying their flight is crucial for translating these biological principles into engineering applications. This proposal outlines a comprehensive computational fluid dynamics (CFD) methodology to simulate the flight of a dragonfly. We detail the governing mathematical equations, the selection of an appropriate numerical scheme, and a step-by-step implementation plan. The proposed framework employs the incompressible Navier-Stokes equations, a pressure-based coupled solver, the Transition SST turbulence model to capture laminar-to-turbulent flow transition, and a dynamic meshing strategy to handle the complex flapping wing kinematics. This research aims to establish a robust and validated simulation model capable of revealing the key fluid dynamics mechanisms, such as leading-edge vortex (LEV) formation and wing-wake interaction, that govern dragonfly flight.
\end{abstract}

\section{Introduction}
Nature has perfected flight over millions of years, and insects, particularly the dragonfly (Odonata), represent a pinnacle of aerodynamic efficiency and performance. With their ability to execute high-speed dashes, intricate aerial maneuvers, and stationary hovering, dragonflies far surpass the capabilities of current MAVs. The key to their performance lies in the complex interaction between their unique corrugated wing structure and the unsteady, three-dimensional flow fields generated by their flapping motion.

The goal of this research is to develop a high-fidelity CFD model to simulate and analyze the aerodynamics of dragonfly flight. By accurately replicating the physical phenomena, we can:
\begin{itemize}
    \item Quantify the aerodynamic forces (lift and drag) generated during various flight modes.
    \item Investigate the formation, evolution, and role of critical flow structures, such as the leading-edge vortex (LEV).
    - Analyze the effect of wing kinematics (flapping, pitching, and deviation angles) on overall performance.
    - Inform the design of next-generation bio-inspired MAVs with enhanced efficiency and maneuverability.
\end{itemize}
This document presents the theoretical foundation and practical implementation plan for achieving these goals.

\section{Mathematical Background: Governing Equations}
The aerodynamics of dragonfly flight, being a low-speed phenomenon (typically 1–10 m/s), is governed by the motion of a viscous, incompressible fluid (air). The fundamental mathematical model for this process is the system of **incompressible Navier-Stokes equations**, which express the conservation of mass and momentum for the fluid.

\subsection{Continuity Equation (Conservation of Mass)}
For an incompressible fluid, where the density \( \rho \) is constant, the conservation of mass is expressed as a statement that the velocity field is divergence-free.
\[
\nabla \cdot \mathbf{u} = 0
\]
Here, \( \mathbf{u} \) represents the velocity vector of the fluid. This equation ensures that the mass of fluid entering any given control volume is equal to the mass exiting it.

\subsection{Momentum Equation}
The momentum equation is essentially Newton's Second Law applied to a fluid element, describing how the velocity of the fluid changes in response to various forces.
\[
\rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} \right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{f}
\]
The terms in this equation are:
\begin{itemize}
    \item \( \rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} \right) \): Represents the inertia of the fluid. The term \( \frac{\partial \mathbf{u}}{\partial t} \) is the local acceleration (unsteadiness), and \( (\mathbf{u} \cdot \nabla)\mathbf{u} \) is the convective acceleration, a non-linear term responsible for complex phenomena like vortex formation.
    \item \( -\nabla p \): The pressure gradient force, which drives the flow from regions of high pressure to low pressure.
    \item \( \mu \nabla^2 \mathbf{u} \): The viscous force, representing the internal friction of the fluid, where \( \mu \) is the dynamic viscosity.
    \item \( \mathbf{f} \): Represents external body forces, such as gravity, which are typically negligible in aerodynamic simulations of this scale.
\end{itemize}
The solution of these coupled, non-linear partial differential equations (PDEs) for the velocity \( \mathbf{u} \) and pressure \( p \) provides a complete description of the flow field around the dragonfly.

\section{Numerical Scheme}
Due to the complex geometry of the dragonfly and the unsteady nature of flapping flight, the Navier-Stokes equations cannot be solved analytically. A robust numerical scheme is required to approximate the solution. The choices for discretization, turbulence modeling, and solver algorithm are critical for accuracy and computational feasibility.

\subsection{Physical Model: Transition SST}
The Reynolds number (Re) for dragonfly flight typically ranges from 100 to 10,000. While the flow may be laminar at the lower end of this spectrum, it is likely to transition to turbulence at higher Re. Therefore, a purely laminar model is insufficient.

We will employ the **Transition SST (Shear Stress Transport) model**. This is a four-equation model that builds upon the standard \(k-\omega\) SST turbulence model by adding two additional transport equations for intermittency and the transition onset momentum thickness Reynolds number. This allows the model to accurately predict the location of laminar-to-turbulent transition on the wing surface, which is a critical physical phenomenon in this flow regime. This model provides an excellent balance between physical accuracy and computational cost, making it ideal for this research.

\subsection{Solver and Discretization}
\begin{itemize}
    \item \textbf{Solver Type:} A **pressure-based coupled solver** will be used. This approach solves the momentum and continuity equations simultaneously, leading to robust and efficient convergence for transient, incompressible flow problems like flapping flight.
    \item \textbf{Spatial Discretization:} To ensure high accuracy, especially in resolving fine vortex structures, the **QUICK (Quadratic Upstream Interpolation for Convective Kinematics)** scheme will be used for the momentum equations. This third-order scheme minimizes numerical diffusion compared to lower-order methods. The gradient will be evaluated using the **Green-Gauss Node-Based** method, which is well-suited for hybrid meshes.
    \item \textbf{Temporal Discretization:} A **First-Order Implicit** scheme will be used for time-stepping. While higher-order schemes exist, this method offers superior stability and robustness, which is crucial when dealing with the large mesh deformations introduced by flapping wings. The time step size, \( \Delta t \), will be chosen to resolve one flapping cycle with 120-160 steps, ensuring that the wing kinematics are captured accurately.
\end{itemize}

\subsection{Handling Wing Motion: Dynamic Meshing}
The large-scale motion of the dragonfly wings is the most significant challenge from a numerical standpoint. A **Dynamic Meshing** strategy is employed to handle this.
\begin{figure}[H]
    \centering
    % \includegraphics[width=0.8\textwidth]{placeholder.png}
    \caption{Conceptual diagram of the hybrid dynamic mesh. A structured C-grid (inner zone) moves rigidly with the wing, while an unstructured mesh (outer zone) deforms and remeshes to absorb the motion.}
    \label{fig:mesh}
\end{figure}
The approach is as follows (see Figure \ref{fig:mesh}):
\begin{enumerate}
    \item \textbf{Structured Inner Zone:} A body-fitted, high-quality C-grid of hexahedral cells will be wrapped tightly around the wing. This zone moves rigidly with the wing's prescribed motion, preserving the mesh quality where aerodynamic gradients are steepest.
    \item \textbf{Unstructured Outer Zone:} The space between the inner zone and the far-field boundaries is filled with an unstructured mesh of tetrahedral cells. This zone absorbs the movement of the inner region.
    \item \textbf{Deformation and Remeshing:} The deformation of the unstructured zone is handled by a **spring-based smoothing** algorithm. To prevent excessive cell distortion, a **local remeshing** method is activated. When the skewness of a cell exceeds a predefined threshold, the solver automatically deletes the poor-quality cells and remeshes that local region, ensuring mesh validity throughout the entire flapping cycle.
\end{enumerate}

\subsection{Implementation Architecture and Programming Structure}
While a production-level simulation would leverage highly optimized software like ANSYS Fluent or OpenFOAM, understanding the underlying programming structure is instructive. A custom solver designed for this problem would be modular, separating the physics, numerical methods, and mesh handling into distinct components.

\begin{itemize}
    \item \textbf{Main Controller/Time-Stepping Loop:} The core of the program is a loop that advances the simulation through time. The high-level logic for each time step is as follows:
    \begin{enumerate}
        \item Get the new wing position and orientation from the pre-defined kinematics function for the current time \(t\).
        \item Call the Dynamic Mesh Module to update the computational grid to conform to the new boundary positions.
        \item Call the Coupled Solver Module to solve the Navier-Stokes equations for the current time step.
        \item Write output data (e.g., flow field, forces) for post-processing.
        \item Advance to the next time step \(t + \Delta t\).
    \end{enumerate}

    \item \textbf{Dynamic Mesh Module:} This module implements the logic for updating the mesh.
    \begin{itemize}
        \item \textit{Rigid Body Motion:} It first applies a solid body rotation and translation to the nodes of the inner structured zone, moving them with the wing.
        \item \textit{Spring-Based Smoothing:} It then iterates through the nodes in the deforming unstructured zone, repositioning them based on the displacement of their neighbors, analogous to a network of springs. The displacement of a node is calculated as a weighted average of its neighbors' displacements.
        \item \textit{Local Remeshing:} After smoothing, the module checks the quality of each cell in the deforming zone (e.g., skewness). If a cell's quality falls below a set threshold, it is marked. The solver then locally deletes the marked cells and the surrounding "cavity" of cells and generates a new, high-quality mesh in that local region.
    \end{itemize}

    \item \textbf{Coupled Solver Module:} This is the computational core. It builds and solves the system of algebraic equations that approximates the PDEs.
    \begin{itemize}
        \item \textit{System Assembly:} The module iterates through every cell in the mesh to assemble a single, large, sparse linear system of the form \( [A][x] = [b] \). The vector \( [x] \) contains the unknown solution variables (e.g., \(u, v, w, p\)) for all cells. The matrix \( [A] \) contains the coefficients that link these variables.
        \item \textit{Discretization:} During assembly, the continuous PDE terms are converted into their algebraic counterparts.
        \begin{itemize}
            \item The \textbf{First-Order Implicit} temporal term is calculated as \( (\phi^{n+1} - \phi^n) / \Delta t \), where \( \phi^{n+1} \) is the unknown value at the new time step and \( \phi^n \) is the known value from the previous step.
            \item The convective term is discretized using the \textbf{QUICK} scheme. To find the value of a variable at a cell face, this scheme uses a quadratic interpolation polynomial fitted through two upstream nodes and one downstream node. This requires a larger data stencil than simpler schemes but provides third-order accuracy and minimizes numerical diffusion, which is critical for preserving sharp vortex structures.
        \end{itemize}
        \item \textit{Linear System Solve:} Once the full matrix \( [A] \) and vector \( [b] \) are assembled, an iterative linear solver (e.g., GMRES, BiCGSTAB) is used to find the solution vector \( [x] \).
        \item \textit{Convergence Loop:} This entire assembly-and-solve process is wrapped in an inner loop that runs for a set number of iterations (e.g., 40-50) or until the solution residuals fall below the specified convergence threshold for that time step.
    \end{itemize}
\end{itemize}

\section{Implementation Process}
The end-to-end process for conducting the CFD simulation follows a structured workflow, from initial geometry preparation to final data analysis.

\subsection{1. Geometry Acquisition and Cleanup}
The process begins with obtaining a high-fidelity 3D model of a dragonfly, likely from a 3D scan (e.g., micro-CT) or a public database. This raw geometry will require cleanup using CAD or mesh editing software (e.g., Blender, MeshLab) to repair defects, fill holes, and simplify overly complex features while preserving key aerodynamic structures like wing venation and corrugations.

\subsubsection{Recommended Approach: Idealized Wing Geometry}
For the initial stages of this research, the most practical and efficient strategy is to use an idealized wing geometry rather than a high-fidelity scan of a real dragonfly. This approach significantly simplifies the workflow and avoids major technical hurdles in mesh generation.

Based on the methodology presented in the literature, we will model the wing as a **zero-thickness, flat quarter-ellipse planform**. This can be easily generated in any standard CAD or mesh pre-processing software (e.g., GAMBIT, as used in the reference studies, or modern equivalents like ANSYS SpaceClaim, Salome, or Blender).

The key advantages of this method are:
\begin{itemize}
    \item \textbf{Simplicity:} It bypasses the need for complex and costly 3D scanning equipment and the subsequent time-consuming cleanup of raw scan data (STL files).
    \item \textbf{Meshing Robustness:} Real dragonfly wings are extremely thin. Attempting to mesh such a thin, solid geometry often results in high-aspect-ratio cells on the wing's edges, which can severely degrade mesh quality, leading to simulation instability and inaccurate results. By modeling the wing as a zero-thickness plate, these problematic "thin face" cells are eliminated entirely, as the wing edge is reduced to a simple line.
\end{itemize}

For this study, we will adopt the dimensions used in the reference paper `Flapping wing design for a dragonfly-like micro air vehicle`:
\begin{itemize}
    \item \textbf{Span:} 7 cm
    \item \textbf{Root Chord:} 4 cm
\end{itemize}
This idealized model provides a well-defined, easily reproducible geometry that still captures the essential planform shape for studying the fundamental aerodynamics of flapping flight.

\subsection{2. Domain and Mesh Generation}
A computational fluid domain (e.g., a large cylinder or cube) will be created around the dragonfly model. The mesh generation will follow the hybrid strategy described in the previous section. A fine boundary layer mesh (prismatic cells) will be generated on the wing and body surfaces to ensure the wall y+ is sufficiently low (\(y^+ \approx 1\)) for the Transition SST model to be effective.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\textwidth]{mesh_concept_diagram.png}
    \caption{Conceptual diagram illustrating the CFD meshing principle. The simulation meshes the fluid domain (air) surrounding the wing boundary, not the solid wing structure itself.}
    \label{fig:mesh_concept}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=\textwidth]{mesh_comparison_visualization.png}
    \caption{A three-way comparison of different mesh visualization techniques: (Left) Raw tetrahedral elements, (Middle) Extracted surface boundaries, and (Right) The fully processed and annotated view. This comparison highlights how different processing steps reveal different aspects of the same mesh data.}
    \label{fig:mesh_comparison}
\end{figure}

\begin{figure}[H]
     \centering
     \begin{subfigure}[b]{0.49\textwidth}
         \centering
         \includegraphics[width=\textwidth]{raw_mesh_structure.png}
         \caption{Raw, unprocessed tetrahedral mesh elements.}
         \label{fig:raw_mesh}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.49\textwidth}
         \centering
         \includegraphics[width=\textwidth]{mesh_structure_visualization.png}
         \caption{Processed mesh with color-coding and slices.}
         \label{fig:processed_mesh}
     \end{subfigure}
     \caption{Visualizations of the computational mesh. The raw view shows a subset of the actual tetrahedral elements, while the processed view uses color-coding and cross-sections to distinguish regions and show internal structure.}
        \label{fig:mesh_visualizations}
\end{figure}

\subsection{3. Solver Setup}
In a CFD software package (e.g., ANSYS Fluent, OpenFOAM), the simulation will be configured with the following detailed settings:
\begin{itemize}
    \item \textbf{Physical Models:} The simulation will be set to solve for incompressible flow using the Transition SST model to handle the laminar-to-turbulent flow transition.
    \item \textbf{Numerical Schemes:} The core numerical setup will use a pressure-based coupled solver. Spatial discretization will be handled by the third-order QUICK scheme for momentum and the Green-Gauss Node-Based method for gradients. Temporal discretization will use the robust First-Order Implicit scheme.
    \item \textbf{Boundary Conditions:} A velocity inlet condition will be set on the upstream boundary of the domain, a pressure outlet on the downstream boundary, and symmetry or slip conditions on the side boundaries. A no-slip wall condition will be applied to the dragonfly's body and wing surfaces.
    \item \textbf{Wing Kinematics:} The complex flapping, pitching, and deviation motions of the wings will be prescribed using User-Defined Functions (UDFs). These functions will define the angular position of the wings as a function of time, typically using sinusoidal equations to mimic the natural motion.
    \item \textbf{Dynamic Mesh Parameters:} The parameters for the dynamic mesh must be set to ensure stability. This includes defining the rigid body motion for the inner zone and setting the parameters for the deforming outer zone. A low spring constant (e.g., 0.001) is recommended for the spring-based smoothing, and appropriate skewness thresholds will be set to trigger local remeshing.
    \item \textbf{Post-processing:} The vast amount of data generated will be analyzed to extract meaningful physical insights. This includes calculating time-averaged and phase-averaged lift and drag coefficients, visualizing pressure and vorticity fields to understand force generation, and tracking the evolution of key flow structures like the LEV and trailing-edge vortices using tools like ParaView or Tecplot.
\end{itemize}

\section{Example Simulation Output}
To illustrate the expected output from the simulation framework, Figure \ref{fig:sim_results} shows two snapshots of the flow field at different points in the flapping cycle. Such visualizations are critical for analyzing the unsteady aerodynamic phenomena that govern dragonfly flight.

\begin{figure}[H]
     \centering
     \begin{subfigure}[b]{0.49\textwidth}
         \centering
         \includegraphics[width=\textwidth]{scene_00020.png}
         \caption{Time Step 20}
         \label{fig:ts20}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.49\textwidth}
         \centering
         \includegraphics[width=\textwidth]{scene_00040.png}
         \caption{Time Step 40}
         \label{fig:ts40}
     \end{subfigure}
     \caption{Instantaneous flow field results at two different time steps. The horizontal slice shows pressure contours, while the vertical slice displays velocity magnitude. These snapshots illustrate the development of complex flow structures around the flapping wing.}
        \label{fig:sim_results}
\end{figure}

\subsection{4. Simulation, Convergence, and Validation}
The transient simulation requires careful management of time-stepping and convergence criteria to ensure a physically accurate and numerically stable result.
\begin{itemize}
    \item \textbf{Time Step Size:} The time step \( \Delta t \) must be small enough to accurately resolve the wing's motion. A common practice is to divide one full flapping cycle into 120 to 160 time steps. For a typical dragonfly flapping at 40 Hz (a period of 0.025s), this corresponds to a \( \Delta t \) of approximately \(1.7 \times 10^{-4}\) s.
    \item \textbf{Convergence Criteria:} Within each time step, the solver performs inner iterations to converge the non-linear equations. A target of 40-50 iterations per time step is typical, aiming to reduce the continuity residuals to approximately \(10^{-4}\) and the momentum residuals to \(10^{-7}\) or lower to ensure the solution is fully converged before advancing to the next time step.
    \item \textbf{Validation:} A mesh independence study will be performed using at least three different mesh resolutions (e.g., coarse, medium, fine) to ensure the results are not dependent on the grid size. The simulation will be run for several flapping cycles to wash out initial transients and achieve a periodic, stable state. Results for lift and drag will be compared with experimental data from the literature where available.
\end{itemize}

\subsection{5. Advanced Steps (Future Work)}
Once the rigid wing model is validated, future work can explore more complex physics, such as:
\begin{itemize}
    \item \textbf{Fluid-Structure Interaction (FSI):} Couple the fluid solver with a structural solver to account for the passive deformation of the flexible dragonfly wings under aerodynamic loads.
    \item \textbf{Higher-Fidelity Turbulence Models:} Employ Detached Eddy Simulation (DES) for a more detailed and accurate resolution of the turbulent wake structures.
\end{itemize}

\end{document} 